Fermion- and spin-counting in strongly correlated systems in and out of thermal 

equilibrium 



Sibylle Braungardt 1 , Mirta Rodriguez 2 , Aditi Sen(De) 3 , Ujjwal Sen 3 , Roy J. Glauber 4 , and Maciej Lewenstein*' 1 
1 ICFO-Institut de Ciencies Fotdniques, Mediterranean Technology Park, 08860 Castelldefels (Barcelona), Spain 
2 Instituto de Estructura de la Materia, CSIC, C/ Serrano 121 28006 Madrid, Spain 
3 Harish- Chandra Research Institute, Chhatnag Road, Jhunsi, Allahabad 211 019, India 
4 Lyman Laboratory, Physics Department, Harvard University, 02138 Cambridge, MA, U.S.A. 
* ICREA-Institucid Catala de Ricerca i Estudis Avancats, 08010 Barcelona, Spain 

Atom counting theory can be used to study the role of thermal noise in quantum phase transitions 
and to monitor the dynamics of a quantum system. We illustrate this for a strongly correlated 
fermionic system, which is equivalent to an anisotropic quantum XY chain in a transverse field, and 
can be realized with cold fermionic atoms in an optical lattice. We analyze the counting statistics 
across the phase diagram in the presence of thermal fluctuations, and during its thermalization 
when the system is coupled to a heat bath. At zero temperature, the quantum phase transition is 
reflected in the cumulants of the counting distribution. We find that the signatures of the crossover 
remain visible at low temperature and are obscured with increasing thermal fluctuations. We find 
that the same quantities may be used to scan the dynamics during the thermalization of the system. 



I. INTRODUCTION 

In the past decade it became clear that the most im- 
portant challenges of physics of ultracold atoms overlap 
essentially with those of condensed matter physics, and 
concern strongly correlated quantum states of many body 
systems. In fact, ultracold fermionic and bosonic atoms 
in optical lattices mimic strongly correlated systems, that 
can be perfectly described by various Hubbard or spin 
models with rich phase diagrams 

Amazingly, atomic physics may address questions con- 
cerning both static and dynamical properties of such sys- 
tems. In the context of statics, the goal is to quantum 
engineer, i.e. to prepare, or reach interesting quantum 
phases or states, and then to detect their properties. 
Many examples of such exotic phases pertaining to quan- 
tum magnetism based on super-exchange interactions are 
now within experimental reach Also, the signatures 

of itinerant ferromagnetism in the absence of the lattice 
structure have been recently reported for a system of spin 
1/2 fermions Q. 

Despite the progress of experimental techniques, the 
preparation and detection of quantum magnetism is al- 
ways obscured by the unavoidable noise and thermal ef- 
fects. These are particularly important in low dimen- 
sional systems, and, in particular, in one dimension (ID) 
where no long range order can exist at T > 0. It is 
therefore highly desirable to design detection methods 
that will allow to observe the signatures of strong corre- 
lations and quantum phase transitions (QFT) at T > 0. 
The first goal of this paper is to demonstrate that atom 
counting may be used to detect signatures of QFTs at 
T > 0. To this aim we analyze a paradigmatic example 
of a strongly correlated system: a system of fermions in 
a ID optical lattice. 

Remarkable long time scales in ultracold atom exper- 
iments allow to monitor the dynamics of the system di- 
rectly. In the context of dynamics, one goal is to observe 



the time evolution of the system under some perturba- 
tion as the system approaches a stationary state. In this 
context various fundamental questions can be addressed. 
For instance, does the system, which can be very well 
regarded as closed, thermalize after initial perturbation 
(sudden quench) [a-flOf? What is the difference between 
thermal and non thermal dynamics? What kinds of in- 
teresting dynamical processes involving a coupling to a 
specially designed heat bath can be realized? Can one re- 
alize state engineering using open system dynamics [HI? 
The second goal of this paper is thus to study atom count- 
ing during dynamic evolution. In particular, we compute 
the atom counting distributions as a function of time 
when the analyzed ID system of fermions approaches the 
quantum Boltzmann-Gibbs thermal equilibrium state at 
certain T > 0. We show how the thermalization pro- 
cess can be monitored by observing the cumulants of the 
counting distribution. In principle, the methods allows 
thus to distinguish thermal dynamics from non-thermal 
one. 

Counting of particles is one of the most important tech- 
niques of characterization of quantum mechanical states 
of many body systems. Photon counting, whose theory 
was developed in the seminal works 12], allows for the 
full characterization of quantum light sources. More re- 
cently, the counting statistics of electrons has been used 
to characterize mesoscopic devices [l3l - fl9| . In both men- 
tioned cases, the particles considered are non-, or practi- 
cally nonintcracting. In this paper, in contrast, we con- 
sider strongly correlated atomic systems [20]. Counting 
statistics of atoms has been suggested as a technique to 
detect and distinguish various quantum phases of spin 
and fermionic systems [2ll - l25j . Atom counting can be 
realized in several manners (for early experiments, see 
[2o|b One method concerns metastable atoms, such as 
metastable Helium [27| - the atoms are released here from 
the trap, so counting is preceded by essentially ballistic 
expansion of the atomic wave functions. With the recent 
development of high-resolution optical imaging systems, 
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single atoms can be detected with near-unit fidelity on 
individual sites of an optical lattice [28, 29]. This makes 
available the counting distributions of atoms in situ in 
the lattice. On the other hand, spin counting techniques 
30] allow for the measurement of the average and fluc- 
tuations of the spin number also in situ in cold atomic 
samples. These techniques can be extended to account 
for spatial resolution [3l| and give access to the Fourier 
components of the spin distribution [32| . With the help 
of superlattice configurations, one may address the atoms 
locally, probing e.g. every second site [3l|. In this work 
we focus on in situ methods, leaving the discussion of the 
interplay of atomic cloud expansion and atom counting 
to a separate publication. 

Despite the fact that in experimental conditions noise 
(thermal or non-thermal) is always present, so far atom 
counting has been mainly considered at zero temperature 
and in the absence of non-thermal noise (2l| - [25j . In par- 
ticular, in Ref. 2 ij . we have used atom counting theory 
to study a system of fermions in a one-dimensional (ID) 
optical lattice. We have shown that the critical behavior 
of the system, and in particular formation of fermionic 
pairs, is reflected in the cumulants of the counting dis- 
tribution. Here, we consider the counting distribution 
of the same fermionic system, but we now take into ac- 
count the effect of thermal noise, both by considering the 
effects of temperature when the system is in an equilib- 
rium state and by that of thermalization when the system 
is attached to a model heat bath. Fermionic pair break- 
ing induced by thermal noise is clearly reflected in the 
counting distribution function. We find also that the sig- 
natures of the crossover between different phases remain 
visible at low temperature, and we show how they fade 
out as the temperature increases. 

The paper is organized as follows. In Sec. [TTl we pro- 
vide a description of the fermion and spin system that 
we consider. In section Mil we review the counting the- 
ory for a fermionic system, and show how the counting 
distribution can be obtained from a simple recursive for- 
mula. Details of how to derive the counting distribution 
in terms of a generating function are shown in the Ap- 
pendix. In section HVt we study the counting statistics of 
the system at thermal equilibrium at non-zero tempera- 
tures. First, in subsection IIV Al we present the counting 
distribution at zero temperature for reference. Then, in 
subsection IIVB1 we analyze how thermal noise affects 
the atom number distributions, especially in the vicinity 
of the quantum phase transition, or better to say cross- 
over. In Sec. |Vl we calculate the atom number distri- 
butions during a model thermalization process, in which 
the system is coupled to a heat bath via the exchange 
of collective quasi-particles. Such couplings, and the re- 
sulting open system dynamics are not strictly speaking 
local. In Sec. IV Bl we analyze, however, the nature of 
these couplings more closely, and show that they can be 
well approximated by a physically feasible model of local 
exchange of atoms between the system and the reservoir. 
We summarize our results in Sec. I VII 



II. FERMI GAS IN A ID OPTICAL LATTICE 

Quantum degenerate fermionic atoms trapped in op- 
tical lattices |33| may become superfluid if there are at- 
tractive interactions between atoms trapped in two differ- 
ent hyperfine states 34]. Attractive fermions form pairs 
analogous to Cooper pairs in superconductors. A one 
component system of fermions trapped in the same hy- 
perfine state may also become superfluid though not in 
s-wave configurations. Such a system, in the ID case, 
can be described by the Hamiltonian (h = 1) 

N 

H = -J ^(cjcj+i + ic]c) +l + h.c. - 2gc]c j + g). 

(1) 

Here, cj denotes the creation of a fermion on site j, N is 
the number of sites, J is the energy associated to fermion 
tunneling to nearest-neighbor lattice sites, g is propor- 
tional to the chemical potential of the system and 7 ac- 
counts for the formation of pairs within consecutive sites. 
A Fourier transform shows that this corresponds to the 
formation and destruction of pairs of opposite momen- 
tum (see 0, [Hj]). A Bogoliubov transformation diago- 
nalizcs the Hamiltonian in Eq. (TTJ , which can be written 
up to a zero energy shift in terms of the quasiparticle 
excitations d k , 

N/2 N/2 

H=J2H k = ^2E k n d k , (2) 
fc=i fc=i 

where 

h d k = d\d k + $_ k d- k (3) 

d k = u k c k - ivkc_ k> d\=u k c\ + iv k c- k , (4) 

4 = ^Ef=iexp(ij$ fe )c], (5) 

u fc = cos^s v k =sin^, (6) 

E k = JyJ (cos$ fc - g) 2 + 7 2 sin 2 $ fe , (7) 

tanfl fc = 7S '^ fc , (8) 

K cos <Pfc — g ' V / 

and <fr k = 2irk/N. In order to recover the Hamilto- 
nian (TT]), for (cos$ fc — g) < the solution of Eq. (j8]) 
is taken from the (^, ^)-branch of the tangent, whereas 
for (cos $fc — g) < it is taken from the (— ■§, ^-branch. 

In the non-interacting case, one can clearly see that 
there are two different regimes. For 7 = the momentum 
space representation of Eq. (fTJ), H k = 2[cos($fc) — g]c\.c k , 
is recovered up to a constant term. For small transverse 
field j < 1, the energy gap of the particles involved is 
positive. For high transverse field g 3> 1, it is negative 
and it vanishes at the critical point g = 1. It can be seen 
from Eq. ((6|) that the coefficients u\ and v\ change their 
roles at the phase transition such that on one side of the 
critical point, the number operator of the quasiparticles 



3 



d\d k corresponds to c\c kl whereas on the other side it 
corresponds to c k (^ k . Finite interactions 7 between the 
fermions lead to the formation of fermionic pairs within 
consecutive sites but the main character of the phase 
transition at g = 1 remains essentially unchanged . 

Quantum phase transitions are only well defined at 
zero temperature. Thermal fluctuations lead to an expo- 
nential decay of the order parameter and only a crossover 
between phases remain. For the system under consider- 
ation the critical point g — 1 at T = extends for finite 
T to a quantum critical crossover region where the en- 
ergy gap is smaller than the thermal fluctuations, i.e. 
\J(l-g)\<k B T^. 

The system considered here (Eq. ([T])) is also interest- 
ing because it is equivalent to the anisotropic quantum 
XY spin model 35[. Using the Jordan- Wigner transfor- 



mation [20, 36], one can transform it into 



N 

- J £ [(1 + 7)^^+1 + (1 - 7)^+1 + gS' s 



(9) 

where S" are the spin 1/2 operators at site j, J is the 
coupling strength, < 7 < 1 is the anisotropy parameter, 
and g is the parameter of the transverse field. The case 
7 = 1 corresponds to the Ising model in a transverse field. 
For 7 = 0, the system corresponds to the isotropic XY- 
model or XX-model. For this value, the Jordan- Wigner 
transformation is ill-defined and one cannot map it to the 
fermionic Hamiltonian Eq.(JTJ. We study the phase tran- 
sition with respect to the parameter g, where the extreme 
cases 3 = and g = 00, correspond to systems with no 
external field and with no interactions, respectively. The 
phase transition between the states with different orien- 
tations of the magnetization takes place at g — 1. For 
small transverse fields g < 1, the ground state has mag- 
netic long-range order and the excitations correspond to 
kinks in domain walls. For high transverse fields g > 1, 
the system is in a quantum paramagnetic state. 



III. FERMION COUNTING STATISTICS 

Before presenting our calculations of the counting dis- 
tribution for a system of fermions at finite temperature, 
we would like to remind ourselves of some basics of pho- 
ton and atom counting statistics. The theoretical analy- 
sis of the counting process of photons registered on a pho- 
todetector was first developed in 37]. In such a process, 
a photon is annihilated and a photoelectron is emitted. 
This photoemission triggers a further ionization process, 
leading to a macroscopic current that is then measured. 
This theoretical framework can be extended for counting 
atoms directly using multichannel plates or in- situ count- 
ing techniques, both for bosons and fermions [38j . In the 
detection process, the particles are absorbed by the de- 
tector. The counting distribution can thus be derived 
from the master equation that describes the interaction 



between the system and the detector with efficiency e 
(see the Appendix). The probability p(m) of counting m 
particles is given by 



p(m) 



(-iy 



dX 1 



-Q 



A=l 



where we have used the generating function 
Q(A) = Tr(p : e~ XI :). 



(10) 



11) 



Assuming that the counting process is much faster than 
the dynamics of the system, the time independent inten- 
sity registered at the detector is I = k^^Li where 
k = 1 — exp(— et) and r denotes the detection exposure 
time. Using the anticommutation relations for fermions 
we obtain 

N/2 

Q(A) = Tr{p\[{l-\Kclc k ){l-\Kcl k c_ k )).(12) 

fe=i 

The dynamics mixes only k and —k fermionic excitations 
dk so that we can separate the density matrix p = TJ fc pf. 
and neglect the terms which do not conserve the number 
of excitations to obtain 



N/2 



(13) 



k=l 



where 

A k = Tr(p k \u 2 k n d k + v 2 k {d k d\ + d_ fc t? fc )] ) 

B k = Tr( Pk \uldld k dlJ^ k + vld^ k dlJ k dl]).{U) 



We use Eq. (|10p to calculate the counting distribution 
from the generating function in Eq. (|13l) and obtain 



p{m) 



(-l) m d m 
m\ d\ 7 



N/2 

■ [ II (1 - XnA k + \ 2 K 2 B k ) 



k=l 



A=l 



n (15) 

Using the generalized Leibniz rule, we derive [24j a recur- 
rence relation to calculate the counting distribution for a 
system with M + 1 pairs of modes from the distribution 
of a system with M pairs of modes 



p(m,M + l) = y2Vip(m-i,M). 



(16) 



Here Vi denotes the probability of detecting i particles 
in the two modes M + 1 and — (M + 1) that is given by 



Pa = 1 - kAm+i + k 2 B m +i, 
Vi = kAm+i - 2k 2 B m +i, 
V2 = I-V0-V1. 



(17) 



Using the recursive relation Eq. (|16p , the counting distri- 
bution for an arbitrarily large system can be calculated 
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from the counting distributions of a two mode system. 
We thus only need to calculate the expressions and 
£>fc in Eq. (fT4")l and use Eqs. (|16H17p to obtain the count- 
ing distributions of the fermionic system Eq.([T]) with an 
arbitrary number of sites. 

As mentioned above, the fermionic operators are re- 
lated to spin operators by the Jordan- Wigner transform. 
The fermion counting distribution is therefore, up to a 
constant, equivalent to the counting distribution of the 
spins in z-direction in the transverse XY-model in Eq. 
© ■ We can thus use the above to calculate the counting 
distributions of the anisotropic XY model in a transverse 
field for a system of any size N . Experimentally, the 
spin number distribution and its fluctuations can be in- 
ferred from the expectation value and fluctuations of the 
polarization of the light that has interacted with a cold 
atomic sample [30| . This spin polarization spectroscopic 
technique can be also extended to account for spatial res- 
olution of the spin distributions [3l| . 

IV. COUNTING STATISTICS IN THE 
PRESENCE OF THERMAL NOISE 

In real counting experiments, there are typically a va- 
riety of noise sources that may affect the system. In this 
section we study the influence of thermal noise on the 
counting distributions of the ID fermi system in Eq. ([TJ. 
We analyze the counting distributions along the crossover 
between the different regions of the phase diagram. We 
first review the results for the zero temperature case and 
then turn our discussion to the case with thermal fluctu- 
ations. 



A. Counting statistics at zero temperature 

At zero temperature, the ground state of the system 
Hamiltonian Eq. (|2|) is the vacuum state of dk excita- 
tions. The expressions A^ and in Eq. (fT4l are thus 
given by 

Ak = 2kv\ 

B k = n 2 v 2 . (18) 

Inserting these into the equations for the two-mode prob- 
abilities Vi in Eq. (TT71) . we obtain the probabilities of 
finding 0, 1 or 2 particles in a system with one pair of 
modes 

p(0, 1) = 1 - 2k{v\ + K 2 vf) 
p(l, 1) = 2kv\ - 2k 2 v\ 

P (2,1) = k 2 v 2 . (19) 

The counting distribution can now be calculated for an 
arbitrary number of modes using the recurrence relation 
Eq. (flU)) . In fig. [1] a , we plot the counting proba- 
bility distribution for the Ising model (7 = 1) at zero 




FIG. 1: a) Counting probability distribution p(m) of finding 
m particles as a function of m for the fermionic system Eq.JT]) 
with 7 = k=1 and g = 0, TV = 1000 at T = 0. b) Mean fh/N 
(blue squares) and variance a 2 /N (red circles) of the counting 
distribution as a function of the transverse field g at T=0. 

temperature with no transverse field (g = 0) and per- 
fect detection efficiency. We consider a system with zero 
excitations and N = 1000 sites. The probability distri- 
bution is centered around a mean value m = 500 = N/2 
particles and its standard deviation a = 50 such that 
a 2 = N/4. As was observed in [24], the pairing that is 
present in the system hamiltonian Eq. ([1]) only allows for 
the detection of pairs of particles and thus leads to a zero 
probability of finding an odd number of particles. In [24j, 
this splitting of the counting distribution between even 
and odd values was shown to disappear for decreasing 
detector efficiency k. In the next section Sec. IIVB1 we 
use this feature of the counting distribution to study the 
influence of thermal fluctuations on the stability of the 
fermion pairs. In fig. [l]b, we plot the mean rh and vari- 
ance a 2 of the counting distribution for different values 
of the transverse field g. The mean number of particles 
increases with increasing transverse field g. The variance 
is constant with g up to the critical point, when it de- 
creases with increasing g. The phase transition at g = 1 
is clearly visible both in the mean and in the variance. In 
[24j , we studied the behavior of the counting distribution 
for different values of the anisotropy parameter 7 and the 
detection efficiency k. We found that the characteristic 
behavior of the mean and variance as shown in fig[T]b) is 
similar when 7 varies from to 1. We further found that 
the phase transition is visible in the means and variances 
even for small detection efficiencies. In the following we 
consider full detection efficiency (k = 1), as the results 
for smaller efficiencies are similar. 



B. Counting statistics at non-zero temperature 

We now turn our discussion to the case of non-zero 
temperature. The effect of thermal fluctuations in the 
system we consider is two-folded. On the one hand, 
thermal fluctuations allow for breaking of superfluid 
fermionic pairs. On the other hand, the quantum phase 
transition reduces to a crossover between different regions 
of the phase diagram. We will show that both effects are 
visible in the counting distribution functions. 

We consider the counting statistics at finite tempera- 
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ture T using the canonical ensemble, p — ie ' , where 
^ = fe^) the Hamiltonian is given by Eq. @ and the 

partition function Z = Tr ^e - ' 3 ^ . The finite tem- 
perature T determines the average number of quasipar- 
ticle excitations dk- In order to calculate the terms Ak 
and Bf, defined in Eq. ((T4)) . we write = -^e~ l3Hk 

where Zk = T5r f e - ^*^ and we take the trace in the 
basis {|00)|01)|10)|11)}. We obtain 

B k = ^(vl + e-^ul) (20) 

Z k = 1 + 2e-? Ek + e~ 2f}Ek (21) 

For a given value of the transverse field g, we fix the 
temperature and obtain the number AT d = J2k=i ^k °f 
fermionic excitations 

Nt = Tr (p fc r$ . (22) 

As explained above, we use A k and B k to obtain the 
recursive formula for the counting distribution. 

Thermal fluctuations induce the breaking of pairs. 
For increasing temperature, the pairing of fermions 
whose binding energy is proportional to 7 in Eq. ([1} is 
suppressed. This is reflected in the counting distribution 
in such a way that the counting probability for odd 
numbers of particles becomes non-zero. To illustrate 
this, we plot in fig. 2 the probability of counting the 
exemplary odd value of m = 499 particles as a function 
of temperature. As temperature increases, the pairs 
are destroyed and we observe a transition from zero 
probability to a finite value. We compare a system with 
small interaction strength 7 = 0.01 (fig. [2] a ) to the 
case of 7 = 1 (fig. [2]b ). In the insets, we compare the 
counting distribution for each system at zero tempera- 
ture and at higher temperatures. We observe that the 
splitting between even and odd particle numbers disap- 
pears as the temperature increases. Note that here we 
consider a perfect detection process. For lower detection 
efficiency, the splitting is not visible, as was shown in 
Ref. [2J|. For small interaction strength 7, the counting 
distribution is narrower, while higher binding energies 
7 imply broader atom number distribution functions. 
Also, observing the scales of temperatures when the 
counting of odd particles become non-zero, one can infer 
that this temperature is proportional to the parameter 7. 

Let us now turn our discussion to the influence of tem- 
perature on the criticality of the system. As was seen 
above for the case of zero temperature, the phase transi- 
tion is visible in the mean and variance of the distribu- 
tion. This behavior is even more evident in the deriva- 
tives of the mean and variance. In fig. [3J we plot the 
derivative of the means and variances with respect to g 



0.2 




k B T/J 



FIG. 2: Probability of counting an odd number (m=499) of 
particles as a function of T for 7 = 0.01 (fig. a) and 7=1 
(fig. b). The insets show the counting distribution for T = 
and k B T/J = 0.01 in fig. a) and for T = and k B T/J = 0.2 
in fig. b). 



at different temperatures T. One can see how the criti- 
cality is blurred when the temperature is of the order of 
the energies of the system /cgT ~ E^. At high tempera- 
ture, the mean and variance become independent of the 
transverse field value g and take a constant value of 0.57V 
and 0.25./V, respectively. 



V. COUNTING STATISTICS DURING 
THERMALIZATION OF A SYSTEM COUPLED 
TO A HEAT BATH 

The long decoherence times of experiments with ul- 
tracold atoms allow to study the real time quantum dy- 
namics of the system. The dynamics of an open system 
coupled to a heat bath has recently aroused much in- 
terest [TlT | as one can use dissipation for quantum state 
engineering. By tuning the properties of the reservoir, 
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FIG. 3: Derivative of the mean m/N (blue solid line) and the 
variance a 2 /N (red dashed line) of the counting distribution 
of the fermionic system Eq. ([l]) with 7 = 1 as a function of the 
transverse field g. a) T = 0, N d = for all g; b) k B T/J = 0.05 
and N d /N ~ at g = 0; c) fc s T/J = 0.3 and N d /N = 0.03 
at fl = 0; d) k B T/J = 1 and N d /N = 0.27 at g = 0. 



thermalization can drive the system to a steady state 
which has the desired properties and can e.g. be used 
to encode quantum information. Here, we consider the 
thermalization of the system hamiltonian Eq. ([lj , when 
it is coupled to a heat bath. We start from the ground 
state at T = and let the system evolve to the thermal 
Boltzmann-Gibbs equilibrium state. In this sense, we 
analyze the counting statistics in a temperature quench. 
Coupling to the heat bath is described by the quantum 
master equation [39| 



dt 



7o 



P(t) = 



1) 



dkp{t)d{ - -d\d k p{t) - -p{t)$ k d k 



d\p{t)d k - ^d k d\p(t) - \^p{t)d k d\ 



(23) 



where 70 is the coupling strength and N k defined in Eq. 
(|22p. accounts for the mean number of fcrmions in the 
mode k at a certain temperature T. This open system 
dynamics assures that the system approaches thermal 
equilibrium towards the Boltzmann-Gibbs state. 

At this point, we would like to clarify an important 
point in relation to particle counting of a dynamical sys- 
tem. The system is governed by two different dynamic 
processes, one is the coupling to the heat bath described 
by Eq. (|2"3")) , the other one is the detection by particle 
counting described by Eq. (|33)l in the Appendix. We 
assume that the coupling of the system to the heat bath 
occurs on a time scale much slower than the counting 
process. The counting is thus performed in a time inter- 



val in which the coupling to the bath does not affect the 
system, so that it can be considered time independent. 
Below we show how the counting statistics change during 
thermalization of the system with the heat bath. How- 
ever, each of the distributions is registered at the detector 
in a time interval in which no change occurs. 



A. Coupling to the excitations 

In order to calculate the counting statistics of the sys- 
tem coupled to a heat bath, we calculate the terms A k 
and B k as given in Eq. (|14[) . which now depend on time. 
From the master equation (|23[) , the time dependent mean 
excitation number is obtained as 



(n d k (t)) 



3 -7o* 



(24) 



We start with the system initially in the vacuum state 
and use 



{d\d k dl k L 



fc (t)> = (dld k ) t (dl k d- k ) 



(25) 



to calculate the time dependent terms A k (t) and B k (t) 
for a system in a heat bath 



= u\N&l - e-^*) + v l(2 - Ng(l - e -»')), 

Bfc(t) _ „2 «(l-e-^)) 2 
— "fc 4 

~ N ^ ~ e ' 70 *) + 7 °' ))2 ). (26) 



In fig. |4l we plot the derivatives of the mean and vari- 
ance with respect to the transverse field g at different 
times t at a fixed coupling rate 70 = 1 and at fixed tem- 
perature of the bath ksT j J = 0.1. At the initial time 
t = 0, the mean and variance correspond to those of the 
zero excitation state, ground state at zero temperature 
(fig. |4] a ) . The phase transition is clearly visible in the 
derivative both of the mean and the variance. Due to the 
coupling of the system and the bath, already for inter- 
mediate times (see fig. [4] b ), the characteristic behavior 
of the mean and variance in the critical region washes 
out. For long coupling time, as shown in fig. U c , the 
behavior is completely determined by the bath. 

In fig. [5j we plot the mean and variance as a function 
of time t for a system coupled to a heat bath at very 
high temperature /cgT/J = 100. Here, the transverse 
field g is fixed. For no transverse field g — (fig. [5] a 
), both the mean and the variance are constant as the 
coupling increases. In the critical point g — 1 (fig. [5] b 
), the variance is constant and the mean decreases as the 
coupling time increases. For high transverse field g = 2 
(fig. c), the mean decreases until reaching the value of 
0.5./V and the variance increases up to the value 0.25^. 



B. Local representation of the coupling 

The master equation Eq. (|2U)) that we use to describe 
thermalization shows two aspects. On the one hand, it 
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is physical to describe the coupling to the bath in terms 
of an exchange of quasiparticles dk, because the Hamil- 
tonian Eq.([T]) conserves the number of quasiparticle ex- 
citations. On the other hand, it may look non-physical 
because the exchange between the system and the bath 
is non-local. The aim of this section is to show that 
the master equation can be rewritten in terms of local 
fermions c; and in p rinciple it could be realized using 
reservoir designs 



.x 10" 



u v 

c)Jt=10 

1 2 

g 

FIG. 4: Thermalization: Derivative with respect to the pa- 
rameter g of the mean (blue solid line) and variance (red 
dashed line) for 7 = 1 for increasing coupling time with 70=1 
and ksT I J — 0.1. Fig. a) shows the initial time when the sys- 
tem is not coupled to the bath. Fig. b) Jt — 1 and c) Jt — 10. 



At high temperatures and in the absence of a trans- 
verse field (g — 0) at any temperature, the number of 
excitations Nf} in the bath is constant with k. In this 
case, the master equation (f2"3)l in terms of the local op- 
erators &i reads 



dt 
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m)pcjc m - 


^uv (J 


m )(p4 d L - 


PClCm))], 



(27) 



where we define the functions 

F„(!-m)4E^ tM 



F tl (/-m) = l 



(l-m) 



F uv (I - to) = — ^2 u k v L < 



i$ fc ((-m) 



(28) 
(29) 

(30) 
(31) 



which depend on the distance I — m between two sites 
I and to and are related to the correlation length of the 
quasiparticles and the pairs. In fig. [6] and fig. [71 we 
study the behavior of the functions F u , F v and F uv as 
the distance between the sites increases. We plot F u , F v 
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FIG. 5: Thermalization: Mean (blue squares) m/N and 
variance a 2 /N (red circles) for increasing coupling time t 
(7 = 7o = 1). a) g = 0, b) g = 1 and c) g = 2. 
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FIG. 6: Thermalization: F u a), F v b) and F uv c) as a func- 
tion of distance between sites (Z — m) for 7 = 1 and g = 
(blue line/circles), g = 1 (green line/squares) and cj = 10 (red 
line /diamands) . 



and i-F™ for different values of g and 7/ J and show that 
the functions F u , F v have their maximum at zero distance 
and decay rapidly as the distance increases. The function 
F uv , which corresponds to the pair correlations, has its 
maximum at the nearest neighbor term I — m = 1. We 
observe that for large transverse field g ;§> 1 , and 7/ J — > 
0, the only non-zero term corresponds to i^(0) = 1. In 
this case, the XY model behaves like a free fermi gas and 
the master equation ([25)) reduces to 



df 



P(t) = 



To 



(AViv + i)£[4 



1.4 1 -4i 



1 



1 



(32) 



Note that for these parameters, the quasiparticles dk — > 
ct. Thus at high T and high transverse field g the bath 
and the system exchange fermionic particles. 

Another interesting limit occurs at any T when g — >• 
and 7/ J = 1. We see in figure [6] that in this case, the 
functions F u , F v are of order 0.5 for the same site and 
F U) F v and F uv are of the order of ±0.25 for neighboring 
sites. The master equation Eq. (l27|) has contributions 
from exchange of on-site fermions and an additional term 
that corresponds to neighboring particles. Also, there 
is exchange of not only on-site particles and holes but 




1 

0.5 






FIG. 7: Thermalization: F u a), F v b) and F uv c) as a func- 
tion of distance between sites Z — m for 7 = 0.01 and g — 
(blue line/circles), g = 1 (green line/squares) and g = 10 (red 
line/diamands) . 
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also fermionic pairs. This is expected as in this regime 
g -C 7, J the pair creation dominates in Hamiltonian Eq. 

©■ 

For low temperatures and at g =/= 0, the number of 
quasiparticles Nj} is not constant with k and the master 
equation cannot be written in the form Eq. (|27|) . How- 
ever, as N% is small for low temperatures, the non-local 
terms are negligible and the equation as a whole remains 
local. 



VI. CONCLUSIONS 

We have studied the effect of temperature in the count- 
ing distribution of a strongly correlated fermionic system 
which can be mapped to the quantum XY spin model 
with a transverse field. Thermal fluctuations induce pair 
breaking in the superfluid fermionic system. We show 
that this is reflected in the number distribution func- 
tion which becomes non-zero for odd number of parti- 
cles for a temperature proportional to the pair formation 
strength. Also, thermal fluctuations reduce the quan- 
tum phase transition into a crossover between different 
regions of the phase diagram. We have found that at 
low temperatures, the mean and variance of the counting 
distribution reflect the critical behavior at the crossover 
between different phase regimes. This effect is obscured 
with increasing temperature and when the temperature 
is comparable to the eigenergies of the system, the critical 
behavior is blurred from the cumulants of the counting 
distribution. 

Furthermore, we have shown that the number distri- 
bution functions can be used to monitor the quantum 
dynamics of the system. We have studied the thermal- 
ization of the system initially at zero temperature when 
it is coupled to a heat bath at finite temperature. This 
process is analogous to a temperature quench. The tem- 
perature determines the number of delocalized excita- 
tions in the system at equilibrium. For high tempera- 
tures and high transverse fields, the exchange of excita- 
tions between system and bath can be mapped into the 
exchange of local fcrmions. For zero transverse field, we 
have shown that the exchange of local excitations cor- 
responds to the exchange of local particles and nearest 
neighbor pairs. We have assumed that the counting pro- 
cess occurs at a different time scale, much faster than the 
exchange of excitations between the system and the bath. 
We have shown that the mean and variance of the count- 
ing distribution can be used to map the thermalization 
process. 

Appendix: Counting of constant fields 

The counting formula in Eq. (1101) has been derived in 
different ways [H, Here, we review a derivation 

(see e.g. [43]) by modeling the absorption of the particles 



at the detector with the master equation 

p — eapa) — —a) tip — —pa) a, (33) 

where and a are the creation and annihilation operator 
of the particle to be counted. Performing a rotation of 
the density matrix, p(t) = e~i taia p(t)e~i ta ^ a , and using 
the relation 

e iA Be - lA — £>_|_ 7[ a B] + f_ [A , [A, B}] + (34) 
we obtain 

p(t) = eae^" 1 pa) e^' = ee~ et apa t . (35) 
This equation can be solved using perturbation theory 

p{t) = p(0) + I ee~ et ' ap{t')a) . 
Jo 

(36) 

Transforming back the rotation we obtain 

p(t) = e~i t&1& (p{0)+ / ee- £t 'ap(0)a t + ...)e-* tata . 
Jo 

(37) 

Using the cyclic properties of the trace, the probability 
Pm(t) of counting to particles can be written as 

( f * dt'ee~ £f ' ) m 
Pm (t) = Tr[p(0)a tm ^ ^_ e -ri«t« a mi (38 ) 

to! 

This is equal to the normally ordered expression 

p m (t) = (:(!- e -*)- (^p e -d-e-)at a ;> (3g) 

to! 

which holds because 

: (a t a) m e-( 1 - e_et ) ata := a) m : e -(i-*- st )tf& ■ a m 

= at" l e - £ * ata a m . (40) 

We can thus use the generating function formalism in 
Eq. (fT2)) with k — f Q dt'ee~ et = 1 — cxp (— er), where r 
is the aperture time of the detector. 
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